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1.0   Introduction 

This  report  concerns  the  use  of  multiquadric  functions  to 
approximate  scattered  data.  Here  we  deal  with  functions  of  two 
independent  variables,  but  the  methods  are  easily  extendible  to 
arbitrary  dimensions,  and  we  expect  that  many  of  the  conclusions 
will  carry  over. 

The  impetus  behind  our  investigation  is  that  of  obtaining 
surface  approximations  that  are  efficient  in  subsequent 
applications.  That  is,  we  consider  it  to  be  acceptable  to  expend 
considerable  computational  resources  to  obtain  the  approximation 
in  a  preprocessing  step.  Once  obtained,  the  approximation  should 
be  able  to  be  evaluated  fairly  efficiently  such  as  when  it  is  to 
be  used  numerous  times  in  an  application  program. 

The  scattered  data  approximation  problem  is  easily  described 
and  occurs  frequently  in  many  branches  of  science.  The  problem 
occurs  in  any  discipline  where  measurements  are  taken  at 
irregularly  spaced  values  of  two  or  more  independent  variables, 
and  is  especially  prevalent  in  environmental  sciences.  We  will 
suppose  that  triples  of  data,  (x^,yj,Zj),  j=l,  ...,  N  are  given, 
assumed  to  be  measurements  (perhaps  with  error)  of  an  underlying 
function  z=f(x,y).  The  function  f  is  to  be  approximated  by  a 
function  F(x,y)  from  the  given  data.  A  recent  survey  of  such 
methods  is  given  in  [FN91]. 

Multiquadric  functions  were  introduced  for  interpolation  of 
scattered  data  by  Hardy  [HA71];  also  see  [HA92]  for  a  historical 
survey  and  many  references.  The  method  is  one  of  a  class  of 
methods  known  now  as  "radial  basis  function  methods"  that 
includes  other  attractive  schemes  such  as  thin  plate  splines 
[HD71,  DU76,  and  others].  The  basic  idea  of  such  methods  is 
quite  simple,  and  we  describe  it  in  some  generality;  for 
purposes  of  being  definite  it  is  pertinent  to  note  that  for  the 
multiquadric  method  the  radial  function  is  h(d)  =  V(d2+r2) .  In 
general,  suppose  a  function  of  one  variable,  h(d),  where  d 
denotes  distance,  is  given. 

For  interpolation  (that  is,  exact  matching  of  the  given 
data),  a  basis  function,  Bj(x,y)  =  h(d^)  is  associated  with  each 


data  point.  Here  dj  =  V( (x-Xj) 2+(y-yj) 2) ,  distance  from  (x,y)  to 
(x-i/Y-i)'  Thus  each  basis  function  is  a  translate  of  the  radial 
function,  h.  The  approximation  is  a  linear  combination  of  the 
basis  functions,  along  with  some  polynomial  terms  that  may  be 
necessary  in  some  cases,  or  may  be  used  to  assure  that  the 
approximation  method  has  polynomial  precision.   Thus, 

N  M 

(1)   F(x,y)  =  £  ajBj(x,y)  +  S  bjqj(x,y) 
j=l  j=l 

where  {q-;}  is  a  set  of  M  polynomials  forming  a  basis  for 
polynomials  of  degree  <m.  The  coefficients  a^  and  b-:  are 
determined  by  the  linear  system  of  equations  prescribing 
interpolation  of  the  data,  and  exactness  for  polynomials  of 
degree   <m: 

N  M 

Z   a^(xi,yi)    +  S   b^(xi,Yi)    =   zL,    i=l,    ...,    N 
j=l  j=l 


(2) 


N 

2   a^qi(x^,y^)    =    0,    i=l,     ...,    M 


D=l 

For  multiquadric  basis  functions,  this  system  of  equations  is 
known  to  have  a  unique  solution  for  distinct  (x-;,y-;)  data  (see, 
for  example,  [MI86]);  while  m  may  be  taken  as  zero  (no 
polynomial  terms) ,  the  theory  indicates  that  a  constant  term 
should  be  included,  and  we  have  done  so  in  all  our  work.  If 
higher  degree  polynomial  precision  is  desired,  inclusion  of  those 
terms  imposes  no  particular  burden. 

While  interpolation  theory  is  important  and  indicates 
something  about  the  suitability  of  the  class  of  functions  for 
approximation  purposes,  our  emphasis  here  is  on  least  squares 
approximation.  This  implies  using  fewer  basis  functions  than 
there  are  data  points.  In  analogy  with  univariate  cubic  splines, 
it  is  convenient  to  refer  to  the  points  are  which  the  radial 
basis  functions  are  centered  as  "knots",  as  was  done  in  [MF92], 
and  we  do  so  here.  If  a  set  of  knot  points,  (uk,vk) ,  k=l,...,K  , 
with  K<N  have  been  specified,  then  the  problem  of  fitting  a 


multiquadric  function  by  least  squares  is  similar  to  that  of 
solving  the  system  of  equations  corresponding  to  those  above  in 
the  least  squares  sense.  We  give  the  details.  Now,  let  B^Cx^y) 
denote    the    radial    basis    function    associated    with    the    point 

(uk'vk)  '    Bk(x'y)="^(  (x"~uk)  2+(y""vk)  2)  *       Ttie    svstem    of    equations, 
specialized  for  our  case,    is  now  of  the   form 

K  M 

2   afc-Bj^x^yi)    +  2   ^q^^y^    =   zif    i=l,     ...,    N 
k=l  j=l 


(3) 


K 

Z   ak   =    0 


k=l 

There  is  a  question  of  how  to  treat  the  last  equation,  which 
guarantees  polynomial  precision.  In  [FC92]  the  corresponding 
constraint  equations  were  imposed  exactly,  rather  than 
approximately  because  of  physical  considerations.  While  there  is 
not  the  corresponding  physical  situation  here,  we  have  also 
imposed  the  last  equation  as  a  constraint.  This  constraint  can 
be  used  to  reduce  the  size  of  the  system  by  solving  for  aK  in 
terms  of  the  other  a^  and  substituting  into  the  first  set  of 
equations. 

If  the  knot  points  are  a  subset  of  the  data  points,  then  the 
same  theory  that  assures  a  unique  solution  of  the  system  for 
interpolation  also  guarantees  a  unique  solution  of  the  least 
squares  problem.  When  the  knot  locations  may  differ  from  the 
data  points,  the  problem  of  whether  the  coefficient  matrix  is  of 
full  rank  or  not  is  unknown  to  us,  although  we  feel  certain  that 
the  matrix  is  of  full  rank  when  the  knot  points  are  distinct,  and 
have  encountered  no  situations  that  indicate  otherwise. 

In  our  implementation  of  the  algorithms  described  in 
subsequent  sections  we  have  used  a  QRP '  decomposition  of  the 
coefficient  matrix  to  solve  the  least  squares  problem.  This 
provides  a  stable  and  efficient  means  for  solution  of  the  problem 
with  an  indication  if  a  matrix  of  less  than  full  rank  is 
encountered . 

In  order  to  test  the  algorithms  we  have  used  a  number  of 


data  sets.  Several  of  these  are  based  on  previously  published 
and  widely  available  (x,y)  data  sets  and  parent  functions.  We 
have  also  used  a  few  less  readily  available  data  sets  that  we  are 
willing  to  share  with  anyone  interested  in  obtaining  them.  Table 
1  gives  a  summary  of  most  of  the  data  set. 


n.m  This  refers  to  point  set  n  and  function  m  from  [FR82], 
for  n=l,  2,  and  3,  and  m=l,  .  ..,  6.  n=l  is  25  points, 
n=2  is  3  3  points,  n=3  is  100  points.  n=4  refers  to  the 
200  point  data  set  used  in  [MF92].  m=l  is  the  humps 
and  dip  function,  m=2  is  the  cliff,  m=3  is  the  saddle, 
m=4  is  the  gentle  hill,  m=5  is  the  steep  hill,  and  m=6 
is  the  sphere.  In  addition,  m=7  refers  to  the  curved 
valley  function  from  [NI78]. 

GT  This  refers  to  the  thinned  glacier  data  consisting  of 
678  points,  with  certain  contour  lines  removed,  from 
[MF92] . 

GL  This  refers  to  the  thinned  glacier  data  consisting  of 
873  points. 

HF  This  is  the  data  set  from  [MF92]  generated  to  be 
approximately  proportional  to  curvature,  consisting  of 
500  points. 

Table  1:   Data  sets  used  extensively  in  tests 

Section  2  deals  with  a  "greedy"  algorithm  for  determining 
the  location  of  a  reasonable  set  of  knots  for  approximation  of 
given  data  by  a  least  squares  multiquadric  function.  Some 
experiences  with  the  method  are  given.  In  Section  3  we  expand 
the  algorithm  to  include  the  knot  locations  and  the  parameter 
value  of  the  method  as  part  of  the  optimization  process.  Some 
results  and  observations  about  the  process  are  made,  with  the 


optimized  value  of  the  parameter  r  being  of  interest.  The 
occurrence  of  near  multiple  knots  is  a  particularly  interesting 
phenomenon.  In  section  4  we  further  extend  the  algorithm  to 
include  variable  parameter  values  at  the  knots.  The  optimized  rk 
values  and  the  near  multiple  knots  are  again  of  special  interest. 
Finally,  in  Section  5  we  discuss  some  ideas  for  further 
exploration  of   least  squares  multiquadric  approximations. 

2.0  An  Adaptive  Method  for  Knot  Selection 

This  section  describes  a  greedy  method  for  the  selection  of 
knot  locations  for  fitting  surfaces  to  scattered  data  using  a 
least  squares  multiquadric  function.  As  noted  in  the 
introduction,  the  use  of  fixed  knots  and  parameter  value  with  the 
multiquadric  function  results  in  a  linear  system  to  be  solved  in 
the  least  squares  sense.  These  are  solved  using  the  QRP ' 
decomposition.  The  algorithm  was  implemented  in  Matlab1,  giving 
access  to  powerful  matrix-vector  notation  that  simplifies  many 
aspects  of  the  implementation.  In  addition,  Matlab  allows  easy 
interactive  intervention  in  the  process,  with  tabular  and 
graphical  information  being  made  available  as  the  computations 
proceeds.  While  an  efficient  implementation  would  also  provide 
for  updating  the  QRP'  decomposition  as  more  knots  are  added,  we 
have  not  done  this  in  the  experimental  program  since  our 
computational   resources  were   sufficient  to  make   it  unnecessary. 

2.1  The  Algorithm 

The  algorithm  proceeds  as  follows,  with  the  necessary  input 
being  obtained  by  interrogation  of  the  user.  The  description 
given  starts  after  all  input  has  been  obtained. 

a)  The  initial  step  is  to  obtain  the  least  squares  fit  by  a 
constant  function,  the  average  of  the  data  values.  The  two 
data  points  having  maximum  positive  and  maximum  negative 
error  are  taken  to  be  the  first  two  knots,  (u1,v1)  and 
(u2,v2).   The  knot  counter  K  is  set  to  2. 

1  MathWorks,  24  Prime  Park  Way,  Natick.  MA  01760 


b)  The  least  squares  multiquadric  fit  with  K  knots  is  obtained. 
The  residuals  are  computed  along  with  their  rms  value,  the 
approximation  is  evaluated  on  a  grid  of  points,  a  smoothness" 
measure  approximately  equal  to  the  value  of  the  thin  plate 
functional  over  the  region  is  computed,  and  if  the  underlying 
function  is  known  the  rms  error  on  the  grid  is  computed.  These 
values  are  then  output  and  a  perspective  plot  of  the 
approximating  surface  is  given. 

c)  The  maximum  absolute  value  of  the  residuals  is  found  and  the 
location  of  this  residual,  subject  to  the  minimum  knot 
separation  value,  is  taken  to  be  the  next  knot  location 
(uK,vK).  At  this  point  the  algorithm  proceeds  to  step  b 
unless  the  maximum  number  of  knot  locations  to  be  computed 
has  been  reached. 

At  the  termination  of  the  program,  the  user  can  restart  the 
process  with  any  of  the  parameters  changed,  with  any  number  of 
knot  point  locations,  up  to  the  total  number  that  have  been 
computed.  Hard  copy  plots  of  the  surfaces  and  tabular  output  can 
be  obtained. 

2.2   Some  Results 

One  of  the  interesting  aspects  of  the  multiquadric  method 
concerns  appropriate  choice  of  the  parameter,  r.  Initial  advice 
was  to  specify  the  value  in  terms  of  approximate  data  point 
separation  [HA71,FR82],  although  even  in  [FR82]  it  was  clear  that 
the  best  value  was  dependent  on  the  ordinate  data  as  well.  More 
recent  work  [TA85,  CF91]  has  shown  this  to  be  the  case  and  an 
algorithm  for  a  "good"  value  was  given  in  [CF91]. 

While  no  algorithm  was  implemented  to  obtain  the  best  r 
value  for  fixed  knots  found  by  the  adaptive  method  above,  the 
flexibility  of  the  implementation  allowed  for  some  interactive 
experimentation  along  these  lines.  In  most  cases  investigated, 
it  was  found  that  the  value  of  r  used  in  the  process  of  selecting 
the  knot  locations  also  was  very  close  to  the  "best"  value  (that 
is  the  one  that  minimized  the  rms  error  of  the  residuals)  for 


that  particular  set  of  knot  locations.  Exceptions  were  when  the 
number  of  knots  was  quite  small  (5*8) , ■ in  which  case  the 
multiquadric  method  shows  a  striking  affinity  for  best  fits  with 
r  very  close  to  zero.  Of  course,  most  surfaces  with  any 
complexity  cannot  be  fit  well  with  so  few  knots.  Apparently  the 
adaptive  knot  selection  process  is  quit  dependent  on  the  r  value 
used,  at  least  enough  to  rule  out  significant  improvement  by 
changing  r  once  the  knots  have  been  selected. 

While  a  reasonable  a  priori  choice  of  the  parameter  r  in 
this  context  can  be  made,  the  value  of  the  best  r  is  still  an 
open  question,  and  is  not  likely  to  be  resolved  anytime  soon.  As 
is  pointed  out  by  [CF91],  the  parameter  can  be  used  somewhat  like 
a  tension  parameter  (small  values  correspond  to  "tighter" 
surfaces) ,  and  consequently  surfaces  that  involve  steep  gradients 
will  be  approximated  with  less  overshoot  by  selecting  a  small 
value  of  r.  The  tension  effects  are  limited  compared  with  the 
results  that  can  be  obtained  using  thin  plate  splines  with 
tension  (see  [FR85]).  Other  factors  enter  into  the  selection  of 
the  r  value,  however,  since  small  values  also  lead  to  rapid 
changes  in  the  gradient  which  may  be  undesirable. 

One  of  the  parameters  in  the  knot  selection  process  is  the 
minimum  separation  between  knots.  It  has  been  found  that  there 
is  often  an  improvement  by  requiring  some  moderate  separation 
between  knots,  for  example  imposing  a  minimum  separation  of  .1 
or  .2  for  2  0  knots  on  the  [0,1]2  for  point  set  3.  This  tends  to 
distribute  knots  more  uniformly  throughout  the  region,  even  when 
there  are  clumps  of  data.  For  comparison  purposes,  the  rms 
errors  (rmse)  at  the  data  points  and  over  a  20x20  grid  were 
computed  and  are  given  in  Table  2  for  several  data  sets.  All  3.m 
examples  were  with  r=0.3  and  20  knot  points,  while  the  HF  data 
set  used  r=0.2  and  50  knot  points. 


data  set  minsep 


3.1 


3.2 


3.3 


3.6 


HF 


nsep 

rmse(data) 

rmse(grid) 

0.0 

0.0215 

0.0231 

0.1 

0.0262 

0.0281 

0.2 

0.0210 

0.0241 

0.0 

0.0122 

0.0132 

0.1 

0.0162 

0.0168 

0.2 

0.0116 

0.0133 

0.0 

0.0020 

0.0023 

0.1 

0.0014 

0.0018 

0.2 

0.0016 

0.0019 

0.0 

0.0033 

0.0038 

0.1 

0.0026 

0.0031 

0.2 

0.0053 

0.0051 

0.0 

0.0031 

0.0040 

0.05 

0.0030 

0.0037 

0.1 

0.0032 

0.0045 

Table  2:   rms  errors  for  various  separation  distances 

In  the  case  given  in  [MF92]  where  the  data  was  specifically 
generated  to  reflect  the  curvature  of  the  underlying  surface,  the 
knots  computed  by  this  algorithm  tend  to  be  gathered  in  regions 
where  the  density  of  data  points  is  greatest.  Figure  1  gives  the 
results  in  one  case.  It  shows  the  data  points  and  the  subset 
selected  as  knot  points  by  the  greedy  algorithm,  along  with  the 
contours  of  the  parent  surface,  in  part  a.  Here  the  minimum 
separation  distance  of  0.05  was  imposed,  resulting  in  a  more 
regular  distribution  than  when  a  zero  separation  distance  is 
imposed.   In  part  b  the  surface  from  which  the  data  was  sampled 


is  shown.  This  function  is  used  in  later  examples  (function  1 
from  the  table) ;  the  viewpoint  is  from  the  right  center  field. 
In  part  c  the  surface  shown  is  that  constructed  by  least  squares 
fit  using  the  knot  points  in  part  a.  Part  d  shows  the  contours 
of  the  approximating  surface.  Part  a  can  be  directly  compared 
with  Figure  3  in  [MF92],  and  it  is  seen  that  the  distribution  is 
different,  and  in  particular  does  not  have  the  nice  spacing  of 
that  in  [MF91].  Qualitatively  the  knot  locations  given  here  do 
reflect  the  density  of  the  data,  however. 

The  greedy  algorithm  given  here  appears  to  be  potentially 
useful  for  many  problems  where  data  subject  to  error  is  available 
and  the  surface  must  be  apprc  iimated  using  an  approximation  that 
is  computationally  as  efficient  as  possible.  A  problem  which  we 
have  considered,  but  which  needs  additional  attention,  is  that  of 
when  enough  knots  have  been  generated  so  that  the  behavior  of  the 
underlying  surface  is  captured  without  undue  influence  by  the 
errors  in  the  data.  For  now,  this  is  mostly  an  unexplored  idea, 
and  we  have  more  to  say  about  it  in  Section  5. 

3.0  Variable  Knot  Locations  and  Multiquadric  Parameter 

While  the  adaptive  method  discussed  in  the  previous  section 
seemed  to  perform  reasonably  well,  it  was  felt  prudent  to  check 
the  performance  of  the  scheme  against  one  which  considered  the 
knot  locations,  along  with  the  parameter  value,  r,  to  be 
variables  over  which  the  minimization  of  rms  errors  at  the  data 
points  was  achieved.  The  function  to  be  minimized  in  this  case 
is  the  same  as  before,  but  here  we  will  state  it  explicitly 
rather  than  in  the  implied  form  where  the  least  squares  solution 
was  that  of  the  overdetermined  system  (3) .  The  minimization 
problem  is 

N       K 
(4)   min  S  [z^   -  2  a^B^^^y^)  -  c]2 
i=l     k=l 

where  the  minimization  is  taken  over  all  (Uj^v^)  ,  r,  the  ak,  and 
c  (with  the  last  equation  of  (3)  imposed  as  a  constraint)  .  As  a 
practical  matter,  for  each  given  knot  configuration  and  r  value, 


the  least  squares  solution  of  (3)  computed  as  a  step  toward  (4)  . 
This  results  in  the  solution  of  a  simpler,  but  equivalent  problem 
since  2K  parameters  are  eliminated  from  (4)  by  imposing  the 
condition  that  the  values  of  the  aj  and  c  be  always  taken  as 
obtained  from  the  least  squares  solution  of  (3).  Hence,  our 
final  process  is  more  properly  written  as 

N       K 
(5)   min  min  2  [z^  -  S  a^Bj^x^y^)  -  c]2 
i=l     k=l 

where  the  inner  minimization  is  over  the  a-;  and  c  (least  squares 
solution  of  (3),  and  the  outer  minimization  is  over  the  knot 
locations  and  the  value  of  the  parameter  r.  The  global  minimum 
of  each  of  the  two  problems  are  clearly  the  same.  Eq.  (5)  is  the 
more  restrictive,  but  any  minimum  of  (4)  is  a  local  minimum  of 
(5) ,  else  a  better  solution  is  attainable  for  (4) .  This  does  not 
imply  that  the  iterative  methods  employed  to  solve  (5)  would  work 
equally  well,  nor  find  the  same  local  minima,  when  applied  to 
(4). 

When  knot  locations  are  allowed  to  differ  from  data 
locations,  the  guarantee  of  full  rank  of  the  coefficient  matrix 
conferred  on  the  system  by  interpolation  theory  no  longer  holds. 
As  we  noted  in  the  Introduction,  this  has  not  posed  any  problems 
in  our  computations. 

3.1  Optimization  Algorithms  and  Initial  Guesses 

We  have  used  two  different  nonlinear  optimization  schemes, 
both  implemented  and  available  as  part  of  Matlab.  One  is  the 
procedure  FMINS  that  is  based  on  a  simplex  procedure  [W085].  The 
other  is  LEASTSQ  that  is  based  on  the  Levinberg-Marquardt 
procedure.  Both  of  the  routines  appear  to  reliably  find  good 
local  minima  that  are  qualitatively  similar,  although  LEASTSQ 
often  finds  a  somewhat  smaller  rms  residual  and  we  have  used  it 
for  most  of  the  results  given  here. 

The  initial  guess  has  a  strong  influence  on  the  solution 
obtained  by  any  nonlinear  optimization  program.  Except  for  a  few 
experiments,  we  have  used  the  results  of  the  greedy  algorithm  in 
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the  previous  section,  with  a  somewhat  judicious  guess  at  the 
value  of  r,  as  the  initial  values  for  the  nonlinear  optimization. 

3.2   Some  Results 

One  of  the  values  of  interest  is  the  optimized  value  of  the 
parameter  r.  For  function  1  the  usual  values  tended  to  be  around 
0.1  to  0.2,  although  in  some  cases  values  outside  that  range  were 
obtained;  the  smallest  rms  errors  were  obtained  in  that  range. 
For  function  2  much  smaller  values  were  obtained,  generally  in 
the  range  less  than  0.05.  For  function  3,  values  in  the  range 
0.20  to  0.30  were  prevalent.  For  function  6  the  value  obtained 
in  the  one  computation  we  carried  out  was  more  than  10.  It  is 
tempting  to  try  to  compare  our  results  with  the  best  values  found 
for  interpolation  by  [CF91],  and  with  their  formula  for 
approximating  the  best  value.  For  the  moment  we  can  say  that  for 
the  most  part  the  data  do  not  seem  contradictory,  although  for 
function  3  our  values  are  somewhat  smaller.  For  function  6,  the 
value  we  obtained  was  in  line  with  computational  experience  in 
[CF91]  in  that  the  value  is  quite  large. 

One  very  interesting  aspect  of  the  results  of  computing 
local  minima  of  (5)  is  that,  with  the  exception  (and  then  not 
always)  of  computations  involving  fewer  than  10  knots,  the 
results  involved  near  repeated  knots,  sometimes  several  different 
pairs  with  20  or  more  knots,  and  sometimes  triples  of  closely 
spaced  knots.  Because  of  the  nonzero  convergence  tolerance  for 
the  optimization  routine,  by  "near  repeated"  knots,  we  mean  those 
that  are  within  a  distance  consistent  with  the  convergence 
tolerance.  In  some  cases  there  were  also  other  knots  within 
distances  of  0.02  or  0.03  for  data  in  the  unit  square. 

The  occurrence  of  near  multiple  knots  suggests  that  the 
method  is  trying  to  adapt  to  some  behavior  of  the  surface  that 
cannot  be  approximated  locally  by  a  single  multiquadric  basis 
function.  The  behavior  of  a  linear  combination  of  multiquadric 
functions  at  points  far  away  is  essentially  the  same  as  a  single 
multiquadric.  Because  of  the  local  extremum  of  the  multiquadric 
function  near  the  knot  point,  it  was  not  immediately  clear  what 
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could  be  achieved  by  a  linear  combination  of  multiquadrics  at 
nearby  knot  points.  Because  of  this,  an  investigation  of  the 
behavior  of  the  surface  defined  by  terms  in  the  approximation  (1) 
corresponding  to  the  near  repeated  knots  was  undertaken,  and  in 
particular,  comparison  with  the  surface  defined  by  a  single 
multiquadric  having  the  knot  point  at  the  average  of  the  repeated 
knots,  with  coefficient  equal  to  the  sum  of  the  coefficients  for 
the  repeated  knots.  Far  away,  the  behavior  of  the  composite 
function  must  be,  and  is  seen  to  be,  essentially  similar  to  a 
single  multiquadric.  In  the  vicinity  of  the  knots  (and  not 
necessarily  just  between  them)  the  behavior  can  be  very 
different. 

Near  multiple  knots  result  in  the  coefficient  matrix  being 
poorly  conditioned,  which  also  allows  for  the  possibility  of 
large  coefficients  in  the  least  squares  solution  of  (3) .  We  are 
unable  to  deduce  for  certain  whether  the  closeness  of  knots  is 
required  in  order  for  the  coefficients  to  become  large,  or 
whether  the  closeness  is  required  to  obtain  the  required  behavior 
in  other  ways.  In  one  case  we  looked  at,  the  knots  are  within 
0.0035  of  each  other,  the  magnitude  of  the  coefficients  is  on  the 
order  of  1125,  and  the  condition  number  of  the  matrix  is  larger 
than  10  ,  some  four  orders  of  magnitude  larger  than  needed  for 
the  magnitude  of  the  coefficients  since  the  data  is  on  the  order 
of  one. 

It  seems  to  be  true  that  the  most  deviant  behavior  of  the 
sum  of  the  near  multiple  terms  occurs  when  the  sum  of  the 
coefficients  for  the  nearby  knots  is  relatively  close  to  zero. 
As  an  example  showing  quite  different  behavior  of  the  sum  of  the 
terms  for  two  nearby  knots  from  that  of  the  average  term,  see 
Figure  2 .  Parts  a  and  b  show  perspective  plots  of  the  two 
surfaces,  while  parts  c  and  d  show  contours  of  the  same  two 
surfaces.  The  deviation  is  striking  and  make  it  seem  reasonable 
that  in  order  to  capture  local  behavior,  multiple  knots  are 
necessary  since  local  behavior  cannot  be  affected  by  basis 
functions  that  are  associate  with  far  away  knots,  and  each  basis 
function  itself  is  locally  a  hyperboloid  (of  one  sheet  -  no 
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saddles)  in  shape. 

Finally,  we  give  the  results  of  optimizing  on  the  knot 
locations  and  the  value  of  r  for  the  50  knot  approximation 
corresponding  to  Figure  1.  The  results  are  shown  Figure  3,  and 
reveal  that  while  many  knots  have  moved  from  their  initial 
positions,  there  density  still  reflects  the  same  general  pattern. 
Noteworthy  is  the  fact  that  there  is  only  one  cluster  of  near 
repeated  knots,  those  being  the  three  at  about  (.44,. 78),  near 
the  dip  in  the  surface.  Those  three  are  clustered  within  a 
distance  of  less  than  0.01,  while  there  is  another  knot  within  a 
distance  of  less  than  0.025.  Finally  we  note  that  the  rms  error 
for  the  surface  was  improved  from  0.0037  to  0.00023  by  the 
optimization  process;   r  changed  from  the  initial  0.2  to  0.2389. 

4.0  Variable  Knot  Locations  Each  with  Variable  Parameter 

The  computational  experience  gained  in  the  variable  knot 
case,  and  especially  the  near  repeated  knot  phenomenon  lead  us  to 
consider  whether  or  not  the  near  multiple  knots  were  occurring 
because  a  single  value  of  the  parameter  at  all  knot  locations  was 
not  necessarily  appropriate.  Thus,  we  modified  the  algorithm  to 
allow  for  an  independent  r  value,  r^,  to  be  associated  with  each 
knot.  The  implications  for  the  rank  of  the  system  is  again  not 
known.  It  is,  however,  easy  to  find  examples  of  different 
parameter  values  that  lead  to  singular  systems  in  the 
interpolation  problem.  We  believe  the  least  squares  problem  is 
more  robust,  and  have  found  no  troublesome  cases  during  our 
investigation . 

4.1  Some  Results 

We  soon  discovered  that  the  use  of  variable  parameter  values 
did  not  alleviate  the  problem  of  near  multiple  knots  in  the 
optimized  approximation.  It  is  interesting  that  when  multiple 
knots  occur,  the  parameter  values  for  those  knots  are  invariably 
very  close  to  having  the  same  value.  In  fact,  our  limited 
experience  seems  to  indicate  that  most  knots  tend  to  have  similar 
values  of  the  parameter,  although  there  are  generally  a  few  that 
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take  on  smaller  values  than  elsewhere. 

Once  again,  the  behavior  of  the  surface  in  the  vicinity  of 
near  multiple  knots  often  reflects  behavior  that  cannot  be  taken 
on  by  a  single  multiquadric.  As  an  example  of  a  different  kind 
of  behavior  than  illustrated  by  Figure  2,  note  that  Figure  4 
shows  another  case  where  the  near  double  knot  results  in  a 
surface  that  resembles  a  quadric  with  a  dimple  in  it.  The  r^ 
values  for  the  two  knots  are  essentially  the  same. 

For  comparison  purposes  between  the  greedy  algorithm,  the 
variable  knot  and  parameter  algorithm,  and  the  variable  knot  each 
with  variable  parameter  value,  we  look  at  a  case  with  a  few 
knots.  In  Figure  5  we  give  the  results  of  the  greedy  algorithm 
for  the  data  set  3-1  with  12  knots  and  an  initial  knot  separation 
of  0.2.  Parts  a-d  are,  respectively,  the  point  and  knot  set,  the 
parent  surface,  the  approximating  surface,  and  contours  of  the 
approximating  surface.  In  Figure  6  the  results  of  the  variable 
knot  and  parameter  algorithm  are  given;  the  initial  values  were 
those  resulting  in  Figure  5.  Parts  a-d  are,  respectively,  the 
point  set,  the  parent  surface,  the  approximating  surface,  and  the 
initial  and  final  knots  locations.  The  improvement  is  clear.  In 
Figure  7  we  see  the  results  of  the  variable  knot,  each  with 
variable  parameter  value.  The  parts  of  the  figure  are  the  same 
as  for  Figure  6.  Here  the  improvement  is  even  more  spectacular. 
The  values  of  the  multiquadric  parameter  and  the  rms  errors  at 
the  data  points  and  on  the  grid  are  given  in  Table  3. 


Algorithm        r  value (s)    rmse(data)     rmse(grid) 

Greedy  0.3  0.0320  0.0341 

Var  kts,  var  r    0.158         0.0101         0.0119 
Var  kts,  var  rk    0.12-0.66     0.0012         0.0023 

Table  3:   Data  set  3-1  with  12  knots,  initial  separation  of  0.2 

The  improvement  in  the  rms  errors  with  variable  knots  is 
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significant  for  this  particular  data  set.  In  situations  where 
the  number  of  knots  is  sufficient  to  give  a  reasonable 
approximation  we  find  the  typical  improvement  in  rms  errors  is 
about  by  a  factor  of  3-10  when  variable  knots  and  parameter 
values  is  allowed,  and  another  factor  of  3-10  when  the  parameter 
values  are  allowed  to  vary.  This  is,  however,  highly  dependent 
on  the  parent  surface,  for  example,  the  cliff  surface 
approximations  are  improved  by  smaller  factors,  while  the  saddle 
surface  approximations  tend  to  be  at  the  upper  end  of  the  scale. 

It  is  interesting  to  compare  the  results  on  this  particular 
example  with  those  of  [CF91]  with  interpolation  to  the  same  data. 
There  it  was  found  that  the  "  est"  value  of  the  parameter  r  (that 
being  about  0.33)  lead  an  approximation  (which  is  the  sum  of  100 
multiquadric  terms)  which  has  an  rms  error  of  0.0026.  It  is 
startling  to  see  that  the  12  term  approximation  derived  using 
variable  knots  each  with  variable  r  has  smaller  rms  error.  We 
have  not  followed  this  line  of  investigation  very  far,  but 
function  2  (cliff)  is  also  approximated  well  using  relatively  few 
knots. 

It  appears  that  the  use  of  variable  knots  can  give  a  greatly 
improved  approximation  when  using  multiquadric  functions  with  a 
fixed  number  of  knots.  When  variable  parameter  values  are 
allowed  the  complexity  of  evaluating  the  approximation  is 
essentially  unchanged  and  seems  to  be  a  worthwhile  improvement 
also.  While  there  is  a  possibility  of  variable  parameter  values 
resulting  in  ill  conditioning  of  the  system,  this  does  not  appear 
to  be  a  real  problem. 

5.0   Conclusions  and  Suggestions  for  Further  Research 

The  methods  we  have  developed  here  appear  to  be  very  useful 
for  the  purposes  we  consider,  that  of  approximating  surfaces  from 
scattered  data  efficiently  for  use  in  subsequent  computations. 
Which  of  the  three  algorithms  one  might  employ  to  obtain  the 
approximation  depends  on  several  matters  that  are  peculiar  to  the 
data  being  approximated,  as  well  as  the  computational 
requirements  and  resources  available.    If  a  reasonable 
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approximation  is  required  with  no  heavy  burden  on  subsequent  use, 
the  greedy  knot  selection  process  will  probably  yield  a  suitable 
result.  If  the  final  use  imposes  a  high  value  on  efficiency  of 
evaluation,  no  doubt  use  of  the  two  schemes  giving  optimized  knot 
locations  will  look  attractive.  From  our  modest  experiments  it 
seems  that  use  of  variable  parameter  values  at  the  knots  adds 
approximation  power  beyond  its  cost. 

To  begin  with,  it  is  desirable  to  carry  out  the 
investigation  with  many  more  sets  of  data.  Exploration  of  the 
process  as  we  have  done  here  is  very  important,  using  known 
underlying  surfaces  being  approximated  for  comparison  purposes. 
However,  ultimately  the  use  of  the  scheme  must  be  for 
approximation  of  data  obtained  experimentally,  or  from 
environmental  measurements.  This  data  is  almost  invariably 
subject  to  error.  While  we  have  does  some  experimentation  with 
such  data  (e,  the  glacier  data) ,  much  remains  to  be  done. 

There  are  a  number  of  directions  in  which  this  research  can 
be  extended.  One  idea  we  have  explored  slightly  is  that  of  using 
some  measure  of  smoothness  of  the  surface,  in  connection  with  the 
rms  error  at  the  data  points,  to  determine  when  to  stop  the  knot 
addition  process  in  the  greedy  scheme.  A  reasonable  stopping 
criterion  is  a  necessity  in  approximating  real-world  data, 
especially  if  the  error  characteristics  are  largely  unknown.  We 
have  computed  the  approximate  value  of  the  thin  plate  functional 
for  these  surfaces  with  the  idea  that  an  significant  increase  in 
the  value  of  the  functional  accompanied  by  only  slight  decrease 
in  the  rms  error  may  signal  that  complexity  is  being  added  to  the 
surface  without  actually  improving  the  fit  to  the  underlying 
surface  by  much.  We  believe  it  will  probably  be  necessary  to 
monitor  the  values  over  some  small  numbers  of  knots,  say  over  5 
or  so  consecutive  numbers  of  knots.  We  intend  to  explore  this  as 
a  potential  stopping  criterion. 

In  certain  sets  of  data  it  may  desirable  to  include  a 
smoothing  term  along  with  the  rms  error  at  the  data  points  as 
part  of  the  objective  function  in  the  knot  location  optimization 
schemes.   One  could  take  the  objective  function  to  be  a  convex 
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combination  of  the  rms  error  and  the  value  of  some  functional 
related  to  smoothness  of  the  surface,  such  as  the  thin  plate 
functional.  There  could  be  several  reasons  for  this  being 
desirable,  but  one  is  that  if  there  are  relative  voids  in  the 
data,  addition  of  a  smoothing  term  would  tend  to  give  some 
control  over  the  behavior  of  the  function  in  such  regions.  There 
are  many  unknown  factors  in  such  a  process.  There  are  numerous 
cases  where  such  objective  functions  have  been  found  useful. 
See,    for   example,     [HS91] . 

The  particular  form  of  the  measure  of  smoothness  probably 
depends  on  the  application,  and  the  use  of  the  thin  plate 
functional,  while  convenient  and  useful  in  many  cases,  may  not  be 
the  proper  one  for  environmental  applications,  for  example.  For 
meteorological  problems  it  has  been  found  that  functionals 
corresponding  to  higher  powers  of  the  Laplacian  seem  to  be 
appropriate  [FR90].  Whatever  the  form  of  the  measure  of 
smoothness,  the  appropriate  choice  of  weighting  between  the  rms 
errors  and  the   smoothness  will   also  have  to  be  discovered. 
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Figure  Captions 


a    b 
All  figures  are  oriented  (sideways)  as 

c    d 

Figure  1:  a)  The  500  data  points  and  the  subset  of  50  of  them 
chosen  as  knot  points  as  generated  by  the  greedy  algorithm  are 
shown.  Minimum  knot  separation  enforced  was  0.05.  Contours  of 
the  parent  function  sampled  for  the  data  are  also  shown.  b)  A 
perspective  plot  of  the  parent  function,  viewed  from  a  point  in 
the  first  quadrant.  c)  A  perspective  plot  of  the  plot  of  the 
approximating  multiquadric  function  computed  as  the  least  squares 
approximation.   d)   Contours  of  the  approximating  function. 

Figure  2 :  The  surface  representing  two  terms  corresponding  to 
two  nearly  repeated  knots  in  a  least  squares  approximation  with 
optimized  knot  locations.  The  region  is  above  the  [0,1]2  square, 
on  which  the  surface  is  sampled.  b)  The  surface  derived  by 
averaging  the  location  of  the  two  knots  and  adding  the 
coefficients.  c)  Contours  of  the  surface  corresponding  to  the 
two  terms  in  part  a.  d)  The  contours  of  the  single  multiquadric 
term  in  part  b. 

Figure  3:  a)  The  parent  surface  sampled  at  the  500  points  show 
in  Figure  la.  b)  The  approximating  least  squares  multiquadric 
with  knots  at  the  initial  guess  points,  as  in  Figure  la.  c)  The 
surface  corresponding  to  the  least  squares  approximation  using 
the  optimized  knot  locations.  d)  The  initial  guesses  and  the 
optimized  knot  locations. 

Figure  4:  The  surface  representing  two  terms  corresponding  to 
two  nearly  repeated  knots  in  a  least  squares  approximation  with 
optimized  knot  locations,  each  with  optimized  multiquadric 
parameter  value;  the  optimized  parameter  values  are  essentially 
the  same.  The  region  is  above  the  [0,1]2  square  on  which  the 
surface  is  sampled.  b)  The  surface  derived  by  averaging  the 
location  of  the  two  knots  and  adding  the  coefficients.  c) 
Contours  of  the  surface  corresponding  to  the  two  terms  in  part  a. 
d)   The  contours  of  the  single  multiquadric  term  in  part  b. 


19 


Figure  5:  a)  The  100  data  points  and  the  subset  of  12  of  them 
chosen  as  knot  points  as  generated  by  the  greedy  algorithm  are 
shown.  Minimum  knot  separation  enforced  was  0.2.  Contours  of 
the  parent  function  sampled  for  the  data  are  also  shown.  b)  A 
perspective  plot  of  the  parent  function,  viewed  from  a  point  in 
the  first  quadrant.  c)  A  perspective  plot  of  the  plot  of  the 
approximating  multiquadric  function  computed  as  the  least  squares 
approximation.   d)   Contours  of  the  approximating  function. 

Figure  6  a)  The  parent  function  which  was  sampled  at  the  100 
points  shown  in  Figure  5a.  b)  The  least  squares  approximation 
constructed  from  the  initial  guess  knot  points,  shown  in  part  d 
as  o's.  c)  The  multiquadric  approximation  constructed  from  the 
optimized  knot  locations  and  (single)  parameter  value.  d)  The 
initial  guess  (o's)  and  optimized  (x's)  knot  locations. 

Figure  7:  a)  The  parent  function  which  was  sampled  at  the  100 
points  shown  in  Figure  5a.  b)  The  least  squares  approximation 
constructed  from  the  initial  guess  knot  points,  shown  in  part  d 
as  o's.  c)  The  multiquadric  approximation  constructed  from  the 
optimized  knot  locations  and  associated  parameter  values.  d) 
The  initial  guess  (o's)  and  optimized  (x's)  knot  locations. 
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